En este notebook nos dedicaremos a estudiar Caminantes Aleatorios. Ellos son la base de movimiento Browniano, tan importante en sistemas físicos y dinámicos.
Ya en el notebook pasado dimos el primer paso a PyCUDA y las grandes ventajas que obtenemos de él. Ya en este notebook podremos hacer algunas cosas interesantes.
La primera pregunta que nos haremos es: ¿cómo hacer que un evento siga una probabilidad $p$, con $0<p<1$?
Pues ya tenemos un generador de números aleatorios. También tenemos condicionales que nos pueden ser útiles. ¿Ya pensaron en la solución? en realidad es muy fácil e intuitiva.
Escribamos un condicional if, de tal forma que al generar un número pseudo-aleatorio $x$, si este es menor a un número $p$, con $0<p<1$, entonces el evento desado pasa.
En otras palabras, si queremos que un evento tenga una probabilidad $0.5$ de pasar, entonces escribiremos un condicional de la forma
x = rand()
p = 0.5
if x <= p:
printf("Hola")
En este ejemplo nuestro evento es imprimir "Hola".
Fácil, ¿no?
[1] Pasemos ahora al primer ejercicio: Escribir una función donde, a partir de una cierta condición inicial, cierta cantidad de threads (a elección del lector) de un paso hacia la izquierda o hacia la derecha ($-1$ o $+1$) con probabilidad $1-p$ y $p$ respectivamente. La función deberá tener como argumento la probabilidad $p$ y deberá regresar el estado resultante de los threads.
In [ ]:
# Primero importamos todos los paquetes necesarios
import pycuda.autoinit
import pycuda.driver as drv
import numpy as np
from pycuda.curandom import rand as curand
from pycuda.compiler import SourceModule
import pycuda.gpuarray as gpuarray
# Primero escribimos el kernel
# El kernel toma como argumentos el arreglo en el cual se lanzarán M números aleatorios con los cuales realizaremos
# los tests probabilisticos (float * proba). También tenemos el arreglo X donde los caminanres se moverán. Finalmente,
# p es la probabilidad de dar un paso a la derecha.
mod = SourceModule("""
__global__ void Caminante(float * proba, float * X, float p){
idx = blockDim.x * blockIdx.x + threadIdx.x ;
if (proba[idx] < p) {
X[idx] += 1 ;
} else {
X[idx] -= 1 ;
}
}
""")
# Define el número de threads M que usaremos y la probabilidad p = 0.5
# Ahora haz dos arreglos en el GPU, uno será la condición inicial de todos los threads.
# El segundo será la probabilidad
# También haz un arreglo semejante al de la condición inicial en el host en el cual copiarás el resultado final
# Obten la función "Caminante" de mod y lanza el kernel
# Copia el resultado al host y finalmente imprimelo
# no te olvides de liberar la memoria
[2] Ahora hagamos una extensión natural de la función pasada.
Tal como en el notebook pasado, ahora completa el código para que el caminante de N pasos hacia la derecha o izquierda. En este ejercicio puedes usar también kernels __device__ para aumentar la eficiencia. También puedes usar la memoria compartida o shared memory.
[3] Finalmente, grafica todas las trayectorias en una gráfica, para distintos valores de $p$ (recomendamos empezar con $p = \frac{1}{2}$). Esto lo puedes hacer con matplotlib tal como vimos en el notebook 01 de la segunda parte.
¿Qué observas?
Regresa a $p = \frac{1}{2}$ e intenta comprender/adivinar/atinarle al comportamiento que sigue la difusión del movimiento de los caminantes aleatorios
Esto lo puedes hacer graficando algunas funciones en la misma gráfica de las distintas trayectorias de los caminantes aleatorios, observando si el comportamiento de las funciones tiene algún parecido con el de los caminantes.
Si ya lo sabías, lo encontraste o ya te cansaste de buscar, te proponemos que veas la imagen siguiente.
Como se puede observar claramente, en este caso la función $\sqrt{n}$ parece describir muy bien el comportamiento de los caminantes aleatorios. Esto no es casualidad y viene de del hecho que el "radio" de un caminante aleatorio se puede relacionar directamente con la deviación estándar que tiene la posición del caminante.
Para una demostración más formal, la página de Mathworld en las referencias contiene una derivación bastante interesante.
Intenta graficar la raíz cuadrada en tu código. ¿no parece estar bien normalizada?
En realidad la función puede aproximarse a $f(n) = A\sqrt{n}$ donde $A$ es un número positivo. Intenta descubir cuál es.
Hint: $A$ está relacionado con el paso que dan los caminantes.
Por ahora hemos dado un primer buen vistazo a los caminantes aleatorios y han aprofundizado sus conocimientos en PyCUDA. Como ya hemos dicho, estos fenómenos podrán ser escalados poco a poco hasta llevar a modelos más complicados, sin perder la esencia básica de los caminantes aleatorios.
Finalmente recomendamos los siguientes ejercicios para que el lector se empape aún más de la magia de los caminantes aleatorios.
[4] Modifica tu caminante para que tenga una probabilidad $r$ de quedarse en el mismo lugar. ¿Cómo cambian los resultados?
[5] Haz lo mismo para caminantes en 2D y caminantes en 3D. [Aquí ya no dibujes las trayectorias como función del tiempo, sino dibuja las trayectorias en el espacio]. Para estos dos casos lo recomendado es usar dos o tres arreglos respecticamente para cada una de las dimensiones (X, Y, Z).
In [ ]: